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Abstract. We calculate global m — 1 modes with low 
pattern speed corresponding to introducing a finite eccen- 
tricity into a protoplanetary disc. We consider disc mod- 
els which are either isolated or contain one or two proto- 
planets orbiting in an inner cavity. Global modes that are 
strongly coupled to inner protoplanets are found to have 
disc orbits which tend to have apsidal lines antialigned 
with respect to those of the inner protoplanets. Other 
modes corresponding to free disc modes may be global 
over a large range of length scales and accordingly be long 
lived. We consider the motion of a protoplanet in the earth 
mass range embedded in an eccentric disc and determine 
the equilibrium orbits which maintain fixed apsidal align- 
ment with respect to the disc gas orbits. Equilibrium ec- 
centricities are found to be comparable or possibly exceed 
the disc eccentricity. We then approximately calculate the 
tidal interaction with the disc in order to estimate the or- 
bital migration rate. Results are found to deviate from the 
case of axisymmetric disc with near circular protoplanet 
orbit once eccentricities of protoplanet and disc orbits be- 
come comparable to the disc aspect ratio in magnitude. 
Aligned protoplanet orbits with very similar eccentricity 
to that of the gas disc are found to undergo litle eccentric- 
ity change while undergoing inward migration in general. 
However, for significantly larger orbital eccentricities, mi- 
gration may be significantly reduced or even reverse from 
inwards to outwards. Thus the existence of global non cir- 
cular motions in discs with radial excursions comparable 
to the semi-thickness may have important consequences 
for the migration and survival of protoplanetary cores in 
the earth mass range. 
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1. Introduction 

The recent discovery of a number of extrasolar giant plan- 
ets orbiting around nearby solar-type stars ( Mayor & 

Send offprint requests to: jcbp@matlis.qmw.ac.uk 



Queloz, 1995, Marcy & Butler 1998, 2000) reveals that 
they have masses that are comparable to that of Jupiter, 
have orbital semi-major axes in the range 

0.04 AU ^ a < 2.5 au, and orbital eccentricities in the 
range 0.0 < e < 0.93. 

Orbital migration originally considered as a mecha- 
nism operating in protoplanetary discs by Goldreich & 
Tremaine (1980) has been suggested as an explanation for 
the existence of giant planets close to their central star 
(eg. Lin, Bodcnheimer, & Richardson 1996). Protoplan- 
etary cores are thought to form at several astronomical 
units and then migrate inwards either before accumulat- 
ing a gaseous envelope and while in the earth mass range 
(type I migration) or in the form of a giant planet (type II 
migration) (see for example Lin & Papaloizou 1993, Ward 
1997). In the former case the interaction with the disc is 
treated in a linearized approximation ( eg. Goldreich & 
Tremaine 1980) while in the latter nonlinearity is impor- 
tant leading to gap formation ( eg. Bryden et al 1999, Kley 
1999). Both types of migration heve been found to occur 
on a timescale more than one order of magnitude shorter 
than that required to form giant planets or the expected 
lifetime of the disc and accordingly the survival of embryo 
protoplanets is in question ( eg. Ward 1997, Nelson et al 
2000 ). Accordingly mechanisms that slow or halt migra- 
tion such as the entry into a magnetospheric cavity close 
to the star ( Lin, Bodenheimer, & Richardson 1996) have 
been suggested. However, the existence of giant planets 
over a range of semi-major axes suggests that a more gen- 
eral mechanism for halting migration should be sought. 

In this paper we shall consider effects produced when 
the large scale motion of the disc gas deviates from that 
of pure circular motion about the central mass as could 
be produced when the disc becomes eccentric as happens 
when it supports a global m — 1 mode. Such modes 
could be produced by disc protoplanet interactions (eg. 
Papaloizou, Nelson & Masset 2001). However, for stan- 
dard parameters, it was found that an instability leading 
to an eccentric disc together with an eccentric protoplanet 
orbit occurred only for large masses exceeding about fif- 
teen Jupiter masses. 
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In this context note these calculations were done for 
laminar model discs with anomalously large viscosity coef- 
ficient and results may be somewhat different for turbulent 
discs. These simulations as well as resonant torque calcu- 
lations appropriate to embedded cores by Papaloizou & 
Larwood (2000) have indicated possible reversal of orbital 
migration associated with eccentric orbits with modest ec- 
centricity comparable to the disc aspect ratio. 

This leads us to study further the dynamics of proto- 
planetary cores, massive enough for tidal interactions with 
the disc to be more important than effects due to gas drag, 
in an eccentric protoplanetary disc. In this paper we con- 
sider large scale m = I modes in gaseous protoplanetary 
discs that correspond to making them eccentric. These 
modes have a low frequency branch for which the disc gas 
follows trajectories differing from Keplerian eccentric or- 
bits by small corrections depending on forces due to disc 
self-gravity and pressure. These modes can be global in 
that they may vary on a length scale comparable to that of 
the whole disc even though it might have a large dynamic 
range. We shall consider a ratio of outer to inner radius 
of one hundred. These modes are of interest because, even 
though no definitive excitation mechanism of general ap- 
plicability has yet been identified, their large scale implies 
a long life time comparable to the viscous time of the disc 
making them of potential interest in Astrophysics ( eg. 
Ogilvie 2001). They have also been considered as a po- 
tential source of angular momentum transport by Lee & 
Goodman (1999) in a tight winding approximation. Fur- 
thermore a disc composed of many stars on near Keplerian 
orbits has been postulated to occur in such objects as the 
nucleus of M31 (Tremaine 1995). 

Here we consider global m — 1 modes for various disc 
models neglecting viscous processes which are presumed 
to act over a longer time scale than that appropriate to 
the phenomena of interest. For global disturbances in pro- 
toplanetary discs of the type we consider, inclusion of self- 
gravity is important. Even though the discs are gravita- 
tionally stable, pressure and self-gravity can be equally 
important on scales comparable to the current radius, r, 
when the Toomre stability parameter Q ~ r/H, with H 
being the disc semi-thickness. This condition is satisfied 
for typical protostellar disc models ( eg. Papaloizou & 
Terquem 1999). In addition to this, disc self-gravity has 
to be considered in the description of the motion of em- 
bedded protoplanets and it generally causes them to move 
in eccentric orbits, with eccentricity comparable to or ex- 
ceeding that in the disc. We find that circularization due 
to tidal interaction with the disc may play only a minor 
role if the local test particle precession frequency relative 
to the disc is large compared to the circularization rate. 
We also find that for modest disc eccentricities comparable 
to the aspect ratio H/r, disc tidal interactions may differ 
significantly from those found in axisymmetric discs. This 
in turn may have important consequences for estimates 



of orbital migration rates for protoplanetary cores in the 
earth mass range. 

In section ^ we give the Basic equations and the lin- 
earized form we use to calculate global m = I modes with 
low pattern speed corresponding to introducing a finite 
disc eccentricity. We go on to describe the disc models 
used which may contain protoplanets orbiting in an in- 
ner cavity. In section |^) we present the results of normal 
mode calculations. We go on to discuss the motion of a 
protoplanet in the earth mass range in an eccentric disc 
in section ^, determining the equilibrium (non precessing) 
orbits which maintain apsidal alignment with the disc gas 
orbits. We then formulate the calculation of the tidal re- 
sponse of an eccentric disc to a low mass protoplanet in 
section ^ determining the time rate of change of the ec- 
centricity and orbital migration rate. We find that aligned 
orbits with very similar eccentricity to that of the gas disc 
may suffer no eccentricity change while undergoing inward 
migration in general. However, when the non precessing 
aligned orbit has a significantly higher eccentricity than 
the disc, as can occur for modes with very small pattern 
speed, orbital migration may be significantly reduced or 
reverse from inwards to outwards for the disc models we 
consider. This finding is supported by a local dynamical 
friction calculation applicable when the protoplanet eccen- 
tricity is much larger than the disc aspect ratio. Finally 
we go on to discuss our results in section ^ 

2. Basic equations and their linearized form for 
large scale global modes with m — 1. 

We work in a non rotating cylindrical coordinate system 
(r, ip) which may initially be considered to be centred on 
the primary star. The basic equations of motion in a two 
dimensional approximation which should be appropriate 
for a large scale description of the disc are taken to be 
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In addition we have the two dimensional form of the con- 
tinuity equation 

d{r^) ^ d{rJ:vr) ^ ^(rSi;^) ^ ^ 
dt dr dip 

Here the velocity is v = (wr, W(^), H = Pdz repre- 
sents a vertically integrated pressure which we assume to 
be a function of the surface density S defined by 

/oo 
pdz. (4) 
-oo 

The sound speed is then given by c = y/dn/dH. The gravi- 
tational potential — ~GM^,/r+^D + ^ext, with G being 



J.C.B. Papaloizou: Protoplanets embedded in eccentric discs 



3 



the gravitational constant, has a point mass contribution 
arising from the central mass Af*, a contribution, due 
to the disc and a contribution, ^ext, due to external pro- 
toplanets. 

The unperturbed disk is axisymmetric with no radial 
motion such that the velocity v — (0, rft) with 17 > 0. In 
equilibrium we then have from equation ([^) 

and the gravitational potential due to the disc is given by 



^D^-G T.K{r,r'ydr', 



with 



Koir,r') 



-\/ (r^ + r'^ — 2rr' cos{(p)) 



dip. 



(6) 



(7) 



We allow for the effects of orbiting external protoplan- 
ets which provide the external potential ^ext ■ In this paper 
we are primarily interested in phenomena which vary on 
a time scale long compared to a local orbital period. Ac- 
cordingly we use the time averaged or secular protoplanet 
perturbing potential which is derived in Appendix 1. For 
a single protoplanet of mass nip and orbital radius Vp, the 
contribution to the external potential is given by 

^ext^-G^Koir,rp). (8) 

In a thin disc of the kind considered here, the contri- 
butions due to pressure and self-gravity are comparable 
and small leading to nearly Keplerian rotation which has 
fl oc r""^/^. When c oc r~^/-^, the disc then has a nearly 
constant putative aspect ratio H/r = c/(rQ), H being the 
putative semi thickness. 

2.1. Linearization 

Here we are interested in global m = 1 modes with 
slowly varying pattern when viewed in the adopted ref- 
erence frame. To study these, we linearize the basic equa- 
tions about the equilibrium state denoting perturbations 
to quantities with a prime. 

As usual we assume that the dependence of all pertur- 
bations (in cylindrical coordinates) on ip and t is through 
a factor exp i (m {(p — ^pt)) . For the slowly varying modes 
we consider Q.p « Q.. 

The linearized forms of equations (|l|) ,(^) and 
then 



dW 



im{n - np)v'^ + —v',. 
and 

im{n - np)T.{w - 0') 



imW 
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(9) 



(10) 



(11) 



Here W = H'/E + $' = S'cVE + 
and k'^ — {2fl/r){d(r'^n)/dr) denotes the square of the 
epicyclic frequency. 

For the gravitational potential perturbation we have 



<I>g^( . Here we shall allow for the contribution 



of the secular effects of other sources such as protoplanets 
to the modes through the gravitational potential pertur- 
bation <i>Lt. 



2.2. Slowly varying modes with m—1 

We consider linear perturbations of an axisymmetric disc 
that correspond to normal modes with m = 1 that make 
it eccentric. The modes we consider are such that flp << 
f2 and the dominant motion is epicyclic. In this case, to 
lowest order, we may neglect the pressure and self- gravity 
term W as well as flp and assume J7 is keplerian, to obtain 
from equation ( |l0|) 



(12) 



Using the above and introducing the radial Lagrangian 
displacement £,r = —iv'^/{Vt — fip), equation (|l]) gives the 
surface density perturbation in the low frequency limit in 
the form 



A-^e{r)) 
dr 



(13) 



with e(r) = ^r/?" being the disc eccentricity. 

The gravitational potential perturbation induced by 
the disc is 



with 



G / Y.'Ki{r,r')r'dr', 



Ki{r,r') 



cos((^) 



(r^ + r'2 — 2rr' cos{(f)) 



dip — 



(14) 



(15) 



where the second term corresponds to the indirect term 
which results from the acceleration of the coordinate sys- 
tem produced by the disc material. 

The existence of a mode with m = 1 in the disc causes 
the orbits of the protoplanets interior to the disc to be- 
come eccentric. Noting the different form of indirect term 
used, the secular perturbing potential derived in Appendix 
1 then gives the contribution of an orbiting protoplanet to 
the external potential perturbation to be 



Gmpirjrp) d 



dr„ 



{K,{r,rpyp) 



(16) 



Here the eccentricity of the companion orbit is related to 
the displacement by Cp = £,rirp)/rp. The total external 
potential perturbation is then found by summing over the 
perturbing protoplanets: 



(17) 
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Eliminating v' from equations (y) and (^^, using 



equation (|l^) for S' and expanding to first order in the 
small frequencies ftp and ujp = CI — k, the fluid element 
orbital precession frequency in the absence of surface den- 
sity or gravitational potential perturbations, one obtains 
a normal mode equation relating e{r),ep and fJp in the 
form (see also Papaloizou, Nelson & Masset, 2001) 



dr 



d (r^^ d[Y.e{r)]\ d{r^^' 



dr 



dr 



-,(18) 



where , c, is the local sound speed. 

In addition we have an equation for each protoplanet 
of the form 



2 {Qp ~ LUp) Q{rp)rpep = - 



d{rH') 
dr 



(19) 



where of course for a particular protoplanet there is no 
self- interaction term in the sum for The inclusion of 
self-gravity in the eigenvalue problem is essential if modes 
with prograde precession frequency are to be obtained. 
For typical protoplanetary disc models, self-gravity can 
be strong enough to induce prograde precession for the 
long wavelength m — 1 modes considered here. 

Normal modes were calculated by discretizing ([l8| ) and 
( |l9| ) and formulating a matrix eigenvalue problem on an 
unequally spaced grid with 200 grid points with intervals 
increasing in geometric progression (see Terqucm & Pa- 
paloizou 2000 for consideration of a related problem) . 



2.3. Disc models 

We present here results obtained for a disc model with 
equation of state given by a polytrope of index n = 1.5. 
The sound speed is given by 



GMJP 



(l - {R,n/r)'°) (l - (r/i?„,0'°) R^n/r.m 



Here the first two factors provide sharp edges at the disc 
inner and outer boundaries and H/r = 0.05 is the disc 
aspect ratio away from these boundaries. 
The surface density is given by 



So(c2)' 



(21) 



The boundary regions are chosen to be of order the 
disc scale height in width and away from these S a r~'^/^ . 
In the above i?i„ is the inner boundary radius which is 
taken to be the unit of length. The unit of mass is the 
central mass and the unit of time is Kll'^ j \fGMl. The 
arbitrary scaling factor Eq was chosen such that in model 
A the total disc mass was 4.0 x 10~^M,. In model B the 
total disc mass was 4.0 x 10"'^. Rout is the outer boundary 
radius, here, in order to study large scale modes, taken as 

100 a„. 



Fig. 1. This figure shows the local test particle orbital 
precession frequency in dimensionless units as a function 
of radius, r in units of i?m in the range 1.1 < r < 100 for 
the two disc models A (lower curve) and B with no interior 
orbiting protoplanets. Apart from small regions near to 
the boundaries, the frequency is negative corresponding 
to retrograde precession. 



The importance of self gravity is measured by the 
Toomre parameter Q — Q.c/ (ttCS). For model A this has a 
minimum value of 5.2 while for the lower mass disc model 
of model B this minimum value is 52. Another quantity 
of interest is the local precession frequency of a test parti- 
cle orbit uopg in the axisymmetric component of the total 
gravitational potential ( see equation ( p9|) below) . This is 
plotted for the two disc models in figure 111. It is small and 
negative in dimensionless units over much of the discs, 
corresponding to retrograde precession, scaling with the 
disc mass as it is determined by the disc self-gravity. As 
we shall see this behaviour provides scope for secular res- 
onance associated with low frequency modes. 



3. Normal modes 

We have calculated the lowest order global normal modes 
for two different configurations involving protoplanets or- 
biting interior to the disc in an inner cavity. The first, 
approximating the Upsilon Andromedae system involves 
two protoplanets with mass ratios 0.00383 and 0.00196 or- 
biting at 0.6i?m and 0.194i?i„ respectively. However, the 
disc normal modes do not change much in character if the 
protoplanet orbital radii are scaled to somewhat smaller 
fractions of i?i„. But the ratio of protoplanet to disc ec- 
centricity in the joint modes changes more significantly. 
The second configuration we consider is a single proto- 
planet with mass ratio 0.002 orbiting at 0.6i?i„. Finally 
we consider a disc with no interior orbiting protoplanets. 
In all of these cases we consider both model A and model 
B discs. 
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Fig. 2. This figure shows the first two modes and the as- 
sociated equilibrium eccentricities as a function of radius 
r in units of Rm for the two protoplanet case with disc 
model A. The two modes (full curves) are normalized so 
that e = 0.1 at r = 1. The lower curve corresponds to 
the mode of highest frequency given in table |l|. The dot- 
ted curves give the associated equilibrium eccentricities 
calculated for a low mass protoplanet orbiting within the 
disc, the lowermost corresponding to the highest frequency 
mode. The dashed curve gives the equilibrium eccentric- 
ity for a low mass protoplanet orbiting within the disc 
calculated for the highest frequency mode but neglecting 
the non axisymmetric component of the disc potential. 
The dot dashed curve gives the corresponding plot for the 
other mode. The latter curves indicate the presence of sec- 
ular resonances at r ^ 1.1 associated with both modes. 

The pattern speeds for the highest frequency normal 
modes and the protoplanet orbital eccentricities occurring 
jointly with the normal modes are given in table 0. For the 
results presented here, the modes are normalized such that 
the disc eccentricity at the inner boundary is 0.1. 

Many of the properties of these modes can be under- 
stood with reference to the local dispersion relation for 
density waves in the low frequency limit (Lin & Shu 1969) 
in the form 

2Vl{uOp-Q.p) = ~2TTGT.\k\+c^k^, (22) 

where k is the radial wavenumber. The above indicates 
that when self-gravity, governed by the first term on the 
right hand side, is unimportant as occurs for either low 
mass discs or large |fc|, fip is negative. On the other hand 
positive Sip can occur for low |fc| or massive discs. Thus in 
that case the lowest order modes (in terms of numbers of 
nodes) should be prograde with more of these existing for 
more massive discs, a trend we find in our results. There 
can also be very low frequency (in magnitude) modes for 
which the self-gravity and pressure terms approximately 
cancel. This occurs when |fc| = 27rGS]/c^ ~ 2/{QH). For 



Fig. 3. As in figure || but for disc model B. In this case 
there is a secular resonance asociated with the highest 
frequency mode at r ~ 1.1. For the other mode this occurs 
at r 1.4. 

Q ^ r/ H as in the models here, this is comparable to the 
radius making a very global mode. 

As we calculate normal modes jointly involving proto- 
planets and disc, the protoplanets have associated eccen- 
tricities. In general the M highest frequency modes pre- 
dominantly involve the protoplanets while the others pre- 
dominantly involve the disc, M being the number of pro- 
toplanets. The two highest frequency modes for the two 
protoplanet case are plotted in figure ^ for disc model A 
and in figure ^ for disc model B. For the second highest 
frequency mode calculated for disc model A, the proto- 
planet eccentricity ratio is 1.38, while for model B it is 
1.88. The latter result is similar to that given by Chiang, 
Tabachnik, & Tremaine (2001) who considered the current 
Upsilon Andromedae system with no disc and concluded 
it was predominantly in this mode. 

The results found here suggest the disc and proto- 
planet orbits were antialigned. With the outermost pro- 
toplanet at O.GRin the disc inner boundary and outer pro- 
toplanet eccentricities are comparable for model A while 
for model B the relative disc eccentricity is 56 percent 
smaller. However, when the outer protoplanet orbits at 
0.5i?i„ the disc inner boundary eccentricity is only 1/3 
that of the protoplanet for model A and 19 percent for 
model B. Apart from the eccentricity scaling relative to 
the interior protoplanets, the spatial form of the eigen- 
function in the disc remains almost identical. 

Apart fom the two modes with highest frequency, other 
normal modes are essentially pure disc modes with only 
small associated protoplanet eccentricities. 

The four highest frequency modes for the one proto- 
planet case with disc model A are plotted in figure ^ and 
in figure |^ for model B. Pure disc modes with no proto- 
planet are plotted in figure ^ for model A and figure |^ for 
model B. The modes develop increasing numbers of nodes 
as their frequency decreases but the mode with lowest fre- 
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Fig. 4. This figure shows the four highest frequency 
modes as a function of radius r in units of i?i„ for the 
one planet model with disc A. The eccentricities are all 
normalized so that e = 0.1 at r = 1. The modes can be 
identified by noting that the number of nodes increases 
with decreasing eigenfrequency. 
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Fig. 5. As in figure ^ but for the one planet model with 
disc model B. 



quency in absolute magnitude can be very global with only 
a few nodes out to r = Rout- We comment that the modes 
are rather non compressive requiring eccentricities of order 
unity to provide Lagrangian changes in surface density of 
order unity. As the motion in the modes in linear theory 
is assumed epicyclic, this is suggestive that the analysis 
should be valid as long as the epicyclic approximation is. 

Having calculated these modes that may be largely 
driven by protoplanets in eccentric orbits or possibly self- 
excited or long lived structures we now go on to consider 
the migration of low mass protoplanets embedded in such 
eccentric discs. 
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Fig. 6. As in figure ^ but for the disc model A with no 
protoplanets. 
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Fig. 7. As in figure ^ but with disc model B 

4. The motion of a low rriEiss protoplanet in an 
eccentric disc 

The evolution of the coplanar orbit of an interior proto- 
planet in the earth mass range embedded in the disc due 
to the gravitational interaction with the slowly precessing 
disc can be found using secular perturbation theory. We 
regard the embedded protoplanet as a test particle evolv- 
ing in a prescribed gravitational potential and neglect the 
influence of the protoplanet on the disc mode. This should 
be reasonable when as here, the angular momentum of the 
protoplanet is significantly less than that of the material 
supporting the normal mode ( eg. Papaloizou et al 2001). 

The orbit evolves under the perturbing potential per 
unit mass, 

V = $D + <^ea:t = $oW + ^-iW C0s(¥. - Qpt). (23) 

Here <l'o(^) is the axisymmetric component of the potential 
but excluding that due to the central mass and $i(r) is 
the radial amplitude of the m ~ 1 component arising from 
the normal mode. 
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Fig. 8. This figure shows the form of the equihbrium 
eccentricity as a function of radius r in units of Rin for 
the modes shown in figure ^ the latter being normahzed so 
that e = 0.1 for r — 1. Curves are associated with normal 
modes according to increasing number of nodes. Thus the 
curve with the most nodes is associated with the normal 
mode with the most nodes. Note the secular resonances 
for r < 5. 




20 ^0 60 80 

Radiu;; 

Fig. 9. As in figure || but for the one planet model with 
disc model B. Equilibrium curves are given for the two 
highest frequency modes only. That associated with the 
highest frequency decays to small values at about r = 5. 
The other curve shows a strong secular resonance at r ~ 
20. 



On performing a time average over the protoplanet orbit 
one obtains the Hamiltonian system: 

dh dTi da dH 

'dt^~'d^' 'dt^'dh' ^ ' 

Here a — w — flpt, with w being the longtitude of pe- 
riapse and Ti. — V — Qph. The specific angular momentum 

is h = J GM^,ap{l — ei), ap is the constant semi-major 
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Table 1. This table gives the protoplanet orbital ec- 
centricities where appropriate together with the pattern 
speeds for the normal modes calculated. The disc model 
used is indicated in the first column, the protoplanet mass 
ratios to the central star in the second and third columns, 
their orbital eccentricities in the fourth and fifth column 
and the mode frequency or pattern speed in the sixth col- 
umn. The modes are normalized such that the disc eccen- 
tricity at the inner boundary is 0.1. A negative protoplanet 
eccentricity then indicates the orbit is antialigned with the 
disc. 
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Fig. 10. As in figure || but for the disc model A with no 
protoplanets. 
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Fig. 11. As in figure p^but for the disc model B with no 
protoplanets. In this case the longest wavelength equilib- 
rium eccentricity indicates antialignment with the disc. 

axis, and Cp is the eccentricity. To leading powers of Cp we 
have 

n = A{ap) + B{ap)el-nph + C{ap)e'^ + epD{ap) cos a.{25) 
Hamilton's equations then give to leading powers of 



da _ 2{B + el{2C - B/2)) 
dt ^GAhap 



_^^_Di^^ (26) 



^jGAUopt 



'dt' 



dn 



Djap) 

yjCM^aptp da y^GM^a 



: sma. 



(27) 



Here the precession rate of the apsidal line induced by 
the axisymmetric potential, to first order in the perturbing 
potential, expressed as an expansion in powers of Cp is 
given by the first term on the right hand of equation (p6|) 
as 



Wo 



{B + el{2C - B/2)) 



(28) 



This may also be written in terms of the axisymmetric 
component of the potential directly as 



2ujpg^/GM*a^ 



(1 



dv? 



4 d4$o 



8a2 du^ ' 



(29) 



where u = 1/r, and r is taken to be equal to ap. 

Equilibrium or steady state solutions of (|2^) and ( p7[ ) 
with Bp and a constant occur when a = or a = tt. Adopt- 
ing the convention that Cp is positive suffices to select one 
of the latter possibilities as (^6|) gives an expression from 
which Cr, can be determined in the form 



D{ap) cos a 



2{B + e2(2C - B/2)) + ^p^GNhap 



(30) 



or equivalently 

ep = D^-P)^^^^ (31) 

il-^pg ~ ^p) ^ GM^Up 

Alternatively one may fix a to be or tt and allow ep to 
change sign as we have done for the normal modes. 

An equilibrium solution so determined corresponds to 
the situation when the protoplanet orbit precesses at the 
same rate, fip, as the mass distribution that produces the 
gravitational potential while maintaining a constant ec- 
centricity. Assuming that ujpg and Op are of comparable 
magnitude, in general for modest eccentricity, 



Cp D{ap)/{ilpy/ GAhttp) = eo 

and is proportional to the magnitude of the nonaxisym- 
metric potential. An exception occurs when 



2B + Vlp^GM^ap = 0. 

In this case, the precession frequency of a free orbit with 
small eccentricity, iOpg, matches that of the nonaxisym- 
metric mass distribution rjp corresponding to a secular 
resonance. When this occurs equation ( |30| ) indicates that 
Cp ~ (eo)^^'^. Accordingly significantly larger equilibrium 
eccentricities are expected close to a secular resonance. We 
comment that for fixed a Cp changes sign as a secular res- 
onance is passed through corresponding to an alignment 
change through a rotation of the axis of the ellipse through 

TT. 

One may also investigate the effect of orbital circular- 
ization by adding a term — ep/|ie| to the right hand side of 
(p7|), where |ie| is the circularization time. In this case one 
can still find an equilibrium but with orbital apsidal line 
rotated. Restricting consideration to the situation away 
from secular resonance, one finds that the equilibrium ec- 
centricity is reduced by a factor yjl -f l/{(ujpg — ^p)'^t1) 
and I sin a I = 1/^1 + (uspg — ^lp)'^t1. Thus when \{u!pg — 
flp)te\ is large, the effect of the circularization term is to 
produce a small rotation of the apsidal line of the orbit. 

As indicated above equilibrium solutions correspond 
to the situation where the apsidal precession of the pro- 
toplanet orbit is locked to that of the underlying nonax- 
isymmetric disc. 

One can find solutions undergoing small librations in 
the neighbourhood of equilibrium solutions (eg. Brouwer 
& Clemence 1961) and when dissipative forces are added 
these may decay making the attainment of equilibrium so- 
lutions natural. As we indicate below, tidal interaction of 
a protoplanet with the disc may produce such dissipative 
forces. Accordingly we shall focus on equilibrium solutions 
in what follows below. 

4-.1. Determination of equilibrium eccentricities 
corresponding to normal modes 

It is a simple matter to determine equilibrium eccentric- 
ities for protoplanets moving under the gravitational po- 
tential appropriate to a normal mode corresponding to 
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an eccentric disc, with pattern speed Qp, of the type cal- 
culated above. We recall the normal mode equation (|l8|) 
which can be regarded as determining the equilibrium disc 
eccentricity in the form 



.(32) 



If the pressure forces are neglected by dropping the 
terms involving c^, this is exactly equivalent to equation 
( ^ ) for determining the protoplanet eccentricity with e(r) 
corresponding to Cp and ujp replaced by cvpg. Thus a local 
protoplanet equilibrium eccentricity appropriate to a disc 
mode or response can be determined by using the disc nor- 
mal mode equation ( |l8|) retaining the term resulting from 
the disc self-gravity but omitting pressure contributions. 

This procedure should be applicable provided the ec- 
centricities are not too large. Equilibrium eccentricities are 
plotted for some of the normal modes we calculated in fig- 
ures H and ^ and also ^ - ^ In figures ^ and || we also plot 
equilibrium eccentricities but calculated with the nonax- 
isymmetric contribution to the gravitational potential due 
to the disc removed which is equivalent to assuming that 
it remains circular. We see that the assumption that disc 
remains circular results in a significantly smaller equilib- 
rium eccentricity for disc model A, particularly for the 
mode with second highest frequency. For the lower mass 
disc model B, differences arising fom assuming the disc re- 
mains circular are less pronounced. We also coment that 
the occurence of secular resonances in the inner parts of 
the disc where r/i?i„ < 10 is common. 

One expects that dissipative torques produced by 
protoplanet-disc tidal interaction may result in the ap- 
proach to such equilibrium solutions from general initial 
conditions. This we now discuss. 



5. Response of an eccentric disc to a lov^r mass 
protoplanet 

To calculate the tidal response of an eccentric disc to 
a perturbing protoplanet it is convenient to introduce a 
new orthogonal coordinate system (a, A), with a being the 
local semi-major axis and A being the orthogonal angu- 
lar coordinate. Then taking pericentre in the disc to lie 
along ip = 0, a — r{l + e cos — e^), with the ec- 
centricity being a function of a. We assume the eccentric- 
ity is small enough that we may work to first order in it 
and also that it has significant changes only on a global 
length scale comparable to r. Then we may make the re- 
placement e(a) = e(r). In this approximation a can be 
expressed in terms of cylindrical coordinates in the form 
a = r(l -I- e cos if) and the orthogonal angular coordinate, 
A is given by 



X = ip - 



/i-dr) sin ip, 
r 



(33) 



where we may leave the integral as indefinite. It is also 
convenient to work in a frame in which the disc flow ap- 
pears stationary. As the pattern speed as measured in an 
inertial frame is very small compared to the disc rotation 
speed, we shall neglect centrifugal and coriolis forces. 

In a general orthogonal coordinate system in which 
the disc appears to be in a steady state, the equations of 
motion for the velocity {va,vx) may then be written ( see 
Appendix 2) 



dt 



iVal' 



d 



d 



^^A|Va||VA|- 



(|V"l) 

dX 



^'iVallVAP- 



IVAI 



2 

iVal dU 



da 



IVal' 



da 



- Va — , 
oa oa 



(34) 



dt 



vx\VX\ 



VAI 



dX 



i;,|Va||VA| 



VAI 



-flVanVAl 



a(|Va|- 
dX 



f |VA|^ 



g(ivAr 

dX 



IvAian a$ 

^9A "'^^'9A' 



(35) 



We perform a response calculation by linearizing about 
a steady state in which the disc gas is taken to be in 
elliptical keplerian orbits with Va = 0, and to first order 
ine,'i;A = (GM, /a) ^/^(l-|-e cos A). In a thin disc forces due 
pressure and self-gravity provide only a small correction 
of order (H/r)'^. We consider the situation when the disc 
eccentricity, e(r), protoplanet eccentricity, ep, and H/r 
can be considered small and comparable in magnitude. 

When considering the response to an embedded pro- 
toplanet, we are interested in responses with a scale H 
and azimuthal mode number m ^ r/H. Here for the time 
being A is regarded as the azimuthal angle. We perform a 
response calculation by linearizing equations (p4) and (35) 



about the steady state described above. We may expand 
the linearized equations in powers of the disc eccentricity. 

Here we shall work only to lowest order and neglect 
terms of order e times smaller than the dominant ones 
which, linearizing (^4|) and (|3^) directly and denoting 
perturbation quantities by a prime, are of order fiw^. In 
this scheme we may replace |Va| and |VA| by unity ex- 
cept where the latter occurs in the combination jVAjwA^, 
with \v\\ here denoting the unperturbed velocity in the 
disc. This is because the contribution of terms of first or- 
der in the eccentricity to this operator leads to quanti- 
ties of order 'me\vx\v'^/r in the linearized equations. For 
m ~ r/H,'^ 1/e these are comparable to the dominant 
terms. All other eccentricity contributions are smaller by 
a factor at most ^ e. 
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We have to first order in eccentricity 
|VA| = (GAf,/a3)i/2 ^i + 2ecos<y9 + J {^dr)cosip^ .(36) 



exp —i{{n — m)Lot — invj + (rti + m)M) 
+ i?2 _ 2rRcos{ip - z/) + &2| 



dMdwdt.{m) 



Thus if we define a new angle M through Af = A — 
(2e + Jifdr) sin A, the operator |VA|wa^, 
becomes = (GM^/a^)^/^^^. The angle M is just the 
mean anomaly as can be seen from the fact that to first 
order in e, 



r = a(l - ecosM), ip = M + 2esinM, 



(37) 



For the motion of a particular disc particle we have M = 
il(t — to), where Iq denotes the time of periastron passage. 

To lowest order in e, the coordinates (a, M) behave just 
like the cylindrical coordinates (r, tp) in the sense that one 
may make the replacements r a, X m the standard 
linearized equations expressed in cylindrical coordinates. 
However, recall that in the forcing potential $p, we must 
make the replacement ( ^t]) which can be thought of eval- 
uation on an eccentric disc orbit. 

5.1. Protoplanet forcing potential 

Neglecting the indirect term, which does not contribute to 
the Fourier components with large m ~ r/H of interest, 
the protoplanet forcing potential is 



$1 



^1^2 _|_ _ 2ri?cos((^ + 6^1 



(38) 



Here, the coordinates of the protoplanet are (i?, v) and h 
is the softening parameter which is introduced to take ac- 
count of the finite vertical thickness of the disc. In general 
h ^ H { eg. Artymowicz 1993). 

For an eccentric protoplanet orbit with semi-major 
axis ap and eccentricity Cp, we have 

R = ap{\ — Cp cosiujt)), ly = Lot + zu + 2ep siii{ujt). (39) 

where the mean motion is lu, we have taken pericentre pas- 
sage to be at i = 0, and w is the longtitude of pericentre 
which is not necessarily zero as we allow for the apsidal 
line of the orbit not to be lined up with that of the disc. 

Using (37) and (|39|), we perform the Fourier decompo- 
sition 



$' = 
p 



E 



Cm,n,ni Gxp i ((n — m)ujt + m{M — vj) + riiM) 



Here, the sum is over positive values of m and both 
positive and negative values of n and ni. The convention 
is that the real part is to be taken to give the physical 
potential. The coefficients Cm,n,ni are real and given by 



The integral is over a 27r cycle in each of the angles 
VJ, M, ut. Note that using ( p7| ) and ( p9[ ) to express the 
disc and protoplanet coordinates makes Cm.n.ni a func- 
tion of a, e, ap, and Cp. 

An important aspect of the linearized problem ex- 
pressed in terms of the coordinates a, M is that the sep- 
arable coordinates are t and M. Accordingly a separable 
response occurs to a combination of terms in the Fourier 
decomposition with fixed (ni +m) = ki and (n — m) — k2. 

Thus for such a particular forcing term with fixed 
(fci,fe) 



*p = fkiM expi{k20Jt + kiM) , 
where 



(41) 



/fci ,fc2 



E 

m,n,ni 



^ni+m,ki ^n~m,k2 ^m,n,nl 

exp(-imti7). (42) 



Here 6 denotes the Kronnccker delta. 

A novel feature is that, unlike in the case of a forced 
axisymmetric disc, a separated forcing potential compo- 
nent in general depends on m. This is because the angle 
between the apsidal lines of the disc and protoplanet or- 
bits can have physical significance. But note that when 
both protoplanet and disc eccentricity are zero only one 
term can survive in (^) when m = ki = — ^2 and then m 
appears only as a redundant complex phase. 



For the general forcing term (41) the pattern speed 



GnipUj 



Q.p = —k2ijj/ki. We shall consider the situation when 
the disc surface density S c>c r~'^/^ as is the case in our 
models away from the boundaries. Then corotation reso- 
nances may be ignored and the main interaction is through 
wave excitation at the Lindblad resonances (eg. Goldre- 
ich & Tremaine 1978). These occur when fci(f2 — 51p) = 
ir^yrn^ with ^ = fcic/(rfi) (eg. Artymowicz 1993). 
For fci > 0, the positive sign applies to the outer Lindblad 
resonance (OLR) and the negative sign to the inner Lind- 
blad resonance (ILR) . In both cases a wave is excited that 
propagates away from the protoplanet. The waves are as- 
sociated with an outward energy and angular momentum 
flux. In the case of the ILR, the background rotates faster 
than the wave so that as it dissipates, energy and angu- 
lar momentum are transferred to the protoplanet orbit. In 
the case of the OLR, the background rotates more slowly 
so that energy and angular momentum are removed from 
'the protoplanet orbit. 

Because the linearized response problem is formally 
identical to the one obtained for an axisymmetric disc 
(a — > r, M ^ (/?) , consistent with he other approxima- 
tions we have made, we can evaluate the outward energy 
flow rate associated with outward propagating waves us- 
ing the approximate expressions developed by Artymow- 
icz (1993) and Ward (1997) which depend only on the 
disc state variables evaluated at the respective resonances. 
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We comment that there are uncertainties associated with 
the use of these expressions in evaluating orbital evolu- 
tion rates especially when they are derived by summing 
torque contributions with varying signs. However, cancel- 
lation effects do not appear to be significant (Artymowicz 
1993, Ward 1997) and as they are supported by general 
considerations, we believe that the main features of the re- 
sults derived to be correct. Following this procedure, the 
outward energy flow rate associated with outward propa- 
gating waves is given by 



dEd _ TT^Sr^ 
where 



(43) 



dr 



However, unlike in an axisymmetric disc, the outward 
angular momentum flow in the disc associated with the 
protoplanet, 
dJd/dt ^ {l/np)dEd/dt. 

One can find the outward angular momentum flow rate 
produced directly by the protoplanet associated with the 
linear response using 



(44) 



the integral being taken over the disc. Recalculation of the 
torque formula then gives 

where 



ki ,k2 



dr 



fir 



with 

9ki,k2 = X! 

m,n,ni 



mSni+mM^n-m,k2Cm,n,nl exp{-imw) .{46) 



The rate of change of protoplanet eccentricity, Cp, induced 
by each term can be found by using the orbit equation 



CriUjJ deri dE dj , 9x_1/9 



(l-e2)3/2 dt 



dt 



(47) 



where E, and J are the orbital energy and angular momen- 
tum of the protoplanet and dJ/dt = ±dJd/dt, dE/dt — 
±dEd/dt^ with the positive sign applying for an ILR and 
the negative sign for an OLR. 

The total rate of change is found by summing contri- 
butions from each resonance occurring for all values of 

(fci, fc2)- 



The eccentricity evolution time-scale, which might be as- 
sociated either with excitation, when it is positive, or 
damping, when it is negative, is then 



(48) 



dcp/dt 



We define the time for the angular momentum to decay 
by a factor of e as 



J 



(49) 



dJ/dr 

such that a positive value means inward migration. 

6. Tidal torques in an eccentric disc 

We have used the above formalism to estimate eccentric- 
ity excitation/damping rates for protoplanets in eccentric 
discs which can be described using the normal modes cal- 
culated above. Before giving details we give a brief sum- 
mary of our results. Although we consider eccentricities 
which may substantially exceed H/r, we shall still sup- 
pose them suflticiently small <~ 0.2 that we may consider 
there to be an equilibrium eccentricity as a function of disc 
radius. Then a torque calculation is characterized by a disc 
eccentricity , e, and an equilibrium eccentricity which we 
may also consider to be the protoplanet eccentricity Cp. 
We shall also for the most part restrict consideration to 
the case when the protoplanet and disc orbits are aligned 
in equilibrium (w = 0.) In general when Cp <^ e we find ex- 
citation of the protoplanet orbit eccentricity, t,, > 0, while 
for Cp ^ e, the eccentricity damps as expected, ie < 0. 
The transition between these regimes occurs when e and Cp 
are approximately equal. For Cp significantly larger than e, 
provided \ {^p — ujpg)te\ is significantly greater than unity, 
there is an equilibrium solution with apsidal line slightly 
rotated from zero. The orbit may then suffer significantly 
reduced or even reversed torques for Cp sufficiently large. 

When Cp is significantly less than a sufficiently large 
disc eccentricity e, the protoplanet orbital eccentricity can 
grow until the orbit ceases to be aligned with the disc and 
precesses through a full {2tt) cycle at which point it then 
damps. Again inward orbital migration may be reduced 
or reversed. Many of these features can be traced to the 
fact that when the equilibrium eccentricity is such that 
Cp — e the situation in many ways replaces the equilibrium 
eccentricity solution Cp — in an axisymmetric disc. This 
we show below 

We use the expressions for r, ip, R, and v, given by 
equations (|37|), ( ^9| ) in the forcing potential (^8|). In the 
case when e = ep, and the protoplanet orbit is almost 
aligned with the disc with very small |ti7|, the strongest 
interaction occurs when r = R, and ip — v. This corre- 
sponds to a — ap, and M = ujt + w. Performing a first 
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Fig. 12. In this figure tj in yr is plotted as a function of 
protoplanet equilibrium eccentricity (negative values cor- 
respond to torque reversal ) for Nt = 1. The short dashed 
curves correspond to an assumed disc eccentricity e — 0.07 
and the short dashed long dashed curves to e = 0.05. The 
disc aspect ratio H/r was taken as 0.07. The crosses were 
obtained for e = 0.05 and H/r — 0.05, while the squares 
correspond to e = 0.07 and H/r = 0.05. In the above cases 
disc and protoplanet orbit apsidal lines were taken to be 
aligned. The diamonds give tj for assumed antialignment 
and e = H/r = 0.07. Note the weakening of the disc in- 
teraction at the higher protoplanet eccentricities due to 
larger relative velocities especially in the antialigned case, 
to the disc. 



Fig. 13. In this figure tf, in yr is plotted as a function of 
protoplanet equilibrium eccentricity (negative values cor- 
respond to circularization) for Nt — 1. The short dashed 
curves correspond to an assumed disc eccentricity e = 0.07 
and the short dashed long dashed curves to e = 0.05. The 
disc aspect ratio H/r was taken as 0.07. The crosses were 
obtained for e = 0.05 and H/r = 0.05, while the squares 
correspond to e = 0.07 and H/r = 0.05. In the above cases 
disc and protoplanet orbit apsidal lines were taken to be 
aligned. The diamonds give for an orbit with apsidal 
line antialigned with that of the disc and e = H/r = 0.07. 
Note the eccentricity grows for protoplanet eccentricity 
less than e in the aligned case while there is always decay 
in the antialigned case. 



order Taylor expansion about the point of maximum in- 
teraction the forcing potential becomes 



p 



(a — flp)^ -|- Aaap sin^(A/ — ut — w) 



62 1 



(50) 



Here |a — Opj and a\M — u}t\ are considered as small com- 
pared to a and terms of order e times these small quanti- 
ties have been neglected. Equation (|5^) is to be considered 
in comparison to the similar expression appropriate to an 
axisymmetric disc 



$' = - 
p 



Grrir, 



(51) 



Given that we have already shown that a, M behave just 
like cylindrical coordinates (r, ip) the tidal torque calcula- 
tion in an eccentric disc with Cp = e should give the same 
results as an axisymmetric disc with Cp = 0. Note too 
that only one term should remain in the sum ([43) such 



that ki = m, and ^2 = —m. Then there is zero rate of 
change of eccentricity. Thus the situation of aligned orbits 
such that Bp — e behaves much like the case with ep = in 
an axisymmetric disc in that it is one of steady eccentric- 
ity. The orbit then migrates at the same rate as a circular 
orbit in an axisymmetric disc. This is essentially what we 
have found on application of the torque formulae. 

6.1. Numerical results 

We have calculated tj/N and te/Nt by summing the con- 
tributions from appropriate resonances. The normalizing 
factor 

Nt = (r/5AU)-i/2(mp/(3Me))-i(I]/56)-i. 

Here the distance is measured in units of 5 AU, the proto- 
planet mass in units of 3 earth masses and the disc surface 
density in units of 56gmcm^^ . We have also taken the cen- 
tral mass to be one solar mass. The softening parameter 
was taken to be & = H/y/2. Our results are consistent with 
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those of Ward(1997) in the hmit where both e and tend 
to zero. We plot tj/Nt as a function of equiUbrium proto- 
planet orbit eccentricity Cp for e = 0.05 and e = 0.07, when 
disc and protoplanet apsidal hnes are ahgned in figure 
Results for aspect ratios 0.05 and 0.07 are presented. 

We plot te for Nt = I in figure |l^. The trends in all 
cases are similar and are that for small Cp the protoplanet 
eccentricity grows in the aligned case while inward migra- 
tion occurs with t j ~ 2 x 10^ yr. The eccentricity growth 
reverses for Cp e indicating an equilibrium in accordance 
with the discussion above. We comment that for very small 
Cp and finite e, dcp/dt oc e. For larger Cp, tj and increase 
and tj eventually changes sign for Cp exceeding ~ 0.1. At 
these eccentricities \te\ ~ 2 x 10'* yr. We make the com- 
ment that the same dimensionless units can be used for 
tj and te as for the disc models introduced in section || 
and thus the same scaling to make results applicable to 
different radii may be used. 

In the antialigned case the disc protoplanet interaction 
is much weaker (see figures |2| and |l^) . Note too that the 
interaction with the disc weakens in general for larger Cp 
because of the larger relative velocity of the protoplanet 
with respect to the disc. This results in larger values of 
\te\ and \tj\. Additional calculations have shown that, as 
expected, these values become independent of the orien- 
tation of the apsidal line when Cp S> e. 

Thus the indications are that for modest eccentricities 
exceeding a few H/r for both disc and protoplanet the 
tidal interaction may differ significantly from the circular 
disc and small protoplanet eccentricity case. Wc now con- 
sider applications to the normal modes calculated in this 
paper. 

The equilibrium protoplanet eccentricities associated 
with the two highest frequency modes calculated in the 
case of disc model A with two interior protoplanets are 
shown in figure ^ while those corresponding to disc model 
B are illustrated in figure ||. These modes are associated 
with significant internal protoplanet eccentricities and can 
be thought of as giving the disc response to external forc- 
ing. Equilibrium eccentricities corresponding to the nor- 
mal modes are also plotted in figures ^ and ^ These are 
generally larger than the disc eccentricities in the inner 
parts of the disc. 

Equilibrium protoplanet eccentricities associated with 
the modes calculated in the case of disc model A with 
one interior protoplanet are shown in figure ^ while those 
corresponding to disc model B are illustrated in figure 
p. Equilibrium eccentricities in the case of isolated disc 
model A are plotted in figure |l^ while those corresponding 
to disc model B are given in figure |ll|. In all of these cases 
the form of the equilibrium eccentricity curve tracks that 
of the corresponding normal mode according to increasing 
number of nodes. Thus the mode with the largest number 
of nodes has associated equilibrium eccentricity with the 
largest number of nodes. 



By comparing the equilibrium eccentricities with their 
corresponding normal modes one sees that the modes with 
the smallest frequencies or pattern speeds in absolute mag- 
nitude tend to have high equilibrium eccentricities several 
times larger than the disc eccentricity. These correspond 
to the longest wavelength curves in figures ^ and ^ with 
corresponding modes plotted in figures ^ and ^. From our 
discussion above these are expected to facilitate high em- 
bedded protoplanet eccentricities. The reason the proto- 
planet eccentricity is significantly larger than the disc ec- 
centricity for these modes is that, for the disc mode the ef- 
fects of the nonaxisymmetric forces due to self-gravity and 
pressure which drive the eccentricity (see equation (p^) 
tend to cancel. However, the protoplanet is subject only to 
self-gravity with no cancelling effects from pressure forces. 
Therefore the equilibrium eccentricity is larger. Note that 
these low frequency modes are essentially disc modes and 
are associated with low interior protoplanet eccentricities 
when the latter are present. Such disc modes may also be 
associated with high embedded protoplanet eccentricities 
through secular resonances ( see section ^ above) . In our 
models such resonances occur when the pattern speed is 
slightly prograde. An example is shown in figure |^. This 
occurs for the one protoplanet model with disc B at about 
20 times the radius of the disc inner edge. Such a reso- 
nance also occurs when disc model A is used but for a 
higher order mode at smaller radii (see figure H) . 

To give numerical examples wc first consider the two 
protoplanet model with disc A as an approximation to the 
Upsilon Andromedae system. Taking the situation rep- 
resented in figure |2| with inner boundary disc eccentric- 
ity 0.1, for a semi-major axis of 2.57 AU for the outer 
protoplanet in the inner cavity, being 0.6 of the inner 
boundary radius, r — 1.5 corresponds to 6.4 AU. For 
the mode with flp — 6.99 x 10^'' and disc model A, 
E = IZQgmcm^^ at r = 1.5. At this radius e = 0.07, 
and the equilibrium eccentricity Bp = 0.1. From figure 
|l^ the orbital migration rate is very small and possi- 
bly outwards. We also find |ie| — lO'*(3Af0/mp) yr and 
\{ijjpg — ilp)ie| ~ 3O.O(M0/mp). Thus in this case the ec- 
centric disc has significant effects on the tidal torques act- 
ing on an embedded protoplanet. 

The eccentricity of the outermost protoplanet orbit in 
the inner cavity would be 0.1. The currently observed 
value which is three or so times larger would be attained 
if the disc inner boundary was scaled to be at twice the 
outer protoplanet semi-major axis. Very similar results for 
the migration and circularization rates to those obtained 
above would still apply. 

As an example to illustrate a case involving a very low 
frequency mode we consider the one protoplanet model 
with disc A. The modes are plotted in figure ^ and the 
equilibrium eccentricities in figure ^. 

For the mode with Vlp — 4.02 x 10^^, we find that 
at r = 10 Bp = 2.5e. Supposing that r — 10 corresponds 
to 5 AU, S = A2(igmcm~^ . From figure |l^ we find that 
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for e = 0.05, there is torque reversal. We also again find 
that \te\ ^ 3 X 10^(3M©/mp) yr and \{ujpg - ftp)te\ ^ 
3Q.0(M©/mp). 

These examples indicate that protoplanets embedded 
in eccentric discs will be in eccentric orbits and that based 
on resonant torque calculations orbital migration slows 
down and may even reverse to become outward when the 
protoplanet eccentricity is sufficiently large. However, we 
should emphasize that these calculations are approximate 
and somewhat uncertain due to the cancellation of torques 
arising at inner and outer Lindblad resonances. To further 
examine the issue of protoplanet disc tidal interaction at a 
large eccentricity ( compared to H/r and e) we present be- 
low a simpler calculation based on local dynamical friction 
which should apply in the appropriate limit and which is 
in essential qualitative agreement with the resonant torque 
calculations. 

6. 2. Dynamical friction calculation 

We here consider the case when Cp significantly exceeds 
H/r but is still also significantly less than unity. This can 
be realized in sufficiently thin discs. In this situation the 
motion of the protoplanet through the disc is supersonic 
so we neglect pressure forces. In addition the scale of the 
response in both space and time becomes local ( even 
though the protoplanet may move globally through the 
disc). For example, for a length scale H, the response time 
scale would be ~ H/{eprVt) ^ . 

Accordingly we work in a reference frame moving in- 
stantaneously with the protoplanet in which the disc ma- 
terial appears to move with velocity v. Adopting local 
Cartesian coordinates, we suppose that the perturbing po- 
tential due to the protoplanet may be written as Fourier 
integral 



$;(r)= j $;(k)exp(ikr)d^k. 



(52) 



Here we shall consider both the two dimensional case 
[N = 2) and the three dimensional case [N — 3). For 
the two dimensional case with $p(r) — —Grrip/y/r'^ -\- 6^, 
$^(k) = -G'mpexp(-&fc)/(27rfc), with k = |k|. In the 
three dimensional case with <i>p(r) = —Gmp/r, ^^(k) = 
-Gmp/(27r2fc2). 

In each case, assuming a local steady state, the velocity 
v' induced by the protoplanet is found from 



vVv' = -V$p. 



We also have 

V • {Dw') = -V(-L>'v), 



(53) 



(54) 



where D denotes the density [N = 3) or the surface den- 
sity {N = 2). In this case where only inertial terms are 
retained the analysis is the same as would be performed 
for coUisionless particles ( eg Tremaine & Weinberg 1984). 



In this context we note that this type of calculation is 
also applicable to the frictional interaction with material 
in planetesimal form as well as with the gas disc provided 
the surface density and disc thickness are appropriately 
specified. 

The rate of change of disc momentum P, which gives 
rise to a frictional force on the protoplanet acting in the 
direction of it's relative velocity, may then be calculated 
from 



P = - / Z^'V*' d r 



(55) 



Performing an integration by parts and working in terms 
of the Fourier transforms of perturbations, one obtains 



P • V -(27r)^ j £)$;(k)ik • v'(k)d^k 

From ( p3| ) one obtains 

-fc2$'(k) 
v'(k).k=— 



(56) 



(57) 



In order to perform the integral ( p6| ) one has to apply a 
Landau prescription by adding an infinitessimally small 
negative imaginary part to the denominator v • k. One 
then finds 



p ^DjGmpf ^^ 



(58) 



where for iV = 3, Q = 4:ln{kmax/kmin), and for N — 2, 
Q = l/{2b). Here {kmax,kmin) are the usual upper and 
lower wavenumber cut offs (eg. Tremaine & Weinberg 
1984). Here reasonable values are k^ax = v^/{Gmp), and 
kmin = The two and three dimensional cases are 

thus of the same form. The logarithmic factor is generally 
of order unity. Comparison of the two and three dimen- 
sional cases suggests that plTi.{kmax / kmin) = S/ {8b). Thus 
the adopted softening parameter b should be somewhat 
smaller than H. 

Using ( pSj ) we may evaluate, remembering that v is the 
relative velocity between disc and protoplanet, the average 
rate of change of angular momentum of the protoplanet 
from 



dJ 



to 
2^ 



7rL>(Gmp)2(r x v) • k 



Qdt, 



(59) 



with the integral being taken round the orbit and k being 
the unit vector in the direction normal to the disc. 

In the two dimensional case for D = T, oc r^^l"^ one 
obtains for small e„ and e = 0. 



d3\ 'KmpUjY.{ap)a-. 
It/ 2be„M? 



(60) 



As this is positive it corresponds to outward migration as 
long as the eccentricity can be maintained. This occurs 
because the disc flow tends to speed up the protoplanet 
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at apocentre where most time is spent. Numerically for 
h = H/^/2 



^ = -2.6 X 10^ 



H 



0.07a, 



/ GpCLp \ 

\ir) 



(61) 



While this agrees in form and is of comparable mag- 
nitude to what is obtained from resonant torque calcula- 
tions (see figure |l^) it is impractical to perform the latter 
at higher eccentricities where better agreement might be 
attained because of the small scale of the interaction (even 
at the eccentricities plotted, over 10^ resonances were in- 
cluded) . 

We may also use the above formalism to calculate the 
mean rate of change of orbital energy for the protoplanet 
to be 



V 



Qdi,(62) 



dE\ _ uj r TTD{Gmp)^v ■ (v - GM^r-^^ 

with E = — '^'^2a"^^ ^'^"^ being the unit vector in the 
azimuthal direction. 

Performing the integration we find for small Cp and as- 
suming E cx r^" that 



1 

'e 



dE 
H 



2ber,M? 



(4.16- 7.17n), 



(63) 



From this we see that, for this circular disc case, the 
mean rate of change of orbital energy increases for n < 0.6. 
This simplified calculation indicates outward migration for 
density profiles that do not increase too rapidly inwards 
in line with the idea that the effect is caused by the inter- 
action at apocentre. 

We further comment that the more sensitive depen- 
dence on the softening parameter means that two dimen- 
sional calculations of the type carried out here, require a 
precise specification of this parameter that correctly repre- 
sents three dimensional effects in order for them to be very 
accurate. Thus two dimensional torque calculations and 
use of torque formulae such as ( |4^ ) suffer from a number 
of uncertainties which can be of comparable importance. 

7. Discussion 

In this paper we have calculated global m — I modes with 
low pattern speed corresponding to introducing a finite 
disc eccentricity. An important aspect was the inclusion of 
self-gravity which is important for the structure of these 
modes as well as determining the motion of embedded 
protoplanets. 

We considered disc models that were isolated or con- 
tained one or two protoplanets orbiting in an inner cavity. 
In all cases global modes were found that could be global 
on scales up to one hundred times the inner cavity ra- 
dius. The modes could be considered as being of a type 
that were strongly coupled to the inner protoplanets or 



essentially free disc modes. In the former case the disc ec- 
centricity could be comparable to that of the protoplanets 
for up to three times the outer protoplanet orbital semi- 
major axis with apsidal line antialigned with that of the 
protoplanet orbits. In the latter case the inner protoplanet 
orbital eccentricities were small compared to that found 
in the disc. 

We went on to discuss the motion of a protoplanet 
embedded in an eccentric disc and determined, initially 
neglecting tidal torques, the equilibrium (non precessing) 
orbits which maintain apsidal alignment with the disc gas 
orbits. Equilibrium eccentricities were found to be compa- 
rable or possibly even exceed the disc eccentricity. In some 
cases secular resonance could occur producing particularly 
large protoplanet eccentricities. 

We then formulated the calculation of the response of 
an eccentric disc to a protoplanet in the earth mass range 
in order to determine the time rate of change of the eccen- 
tricity and orbital migration rate. We found that equilib- 
rium aligned orbits with very similar eccentricity to that 
of the gas disc may suffer no eccentricity change while 
undergoing inward migration in general. This was found 
from the resonant torque calculations but is also expected 
from direct consideration of the equations governing the 
tidal response. However, when the non precessing equilib- 
rium aligned orbit has a significantly higher eccentricity 
than the disc, as can occur generally, but in particular for 
modes with very small pattern speed, orbital migration 
may be significantly reduced or reverse from inwards to 
outwards for the disc models we considered. 

Attainment of high eccentricities in this way typi- 
cally requires the characteristic test particle orbit preces- 
sion frequency or mode pattern speed to significantly ex- 
ceed the characteristic orbital circularization rate, a sit- 
uation more likely for lower mass protoplanets. When 
tidal circularization dominates, the protoplanet equilib- 
rium eccentricity is reduced while the apsidal line be- 
comes significantly inclined to that of the disc. However, 
high protoplanet eccentricities could be excited by grav- 
itational interactions between them (Papaloizou & Lar- 
wood 2000) and under favourable conditions this effect 
could act to counter tidal circularization generating sig- 
nificantly higher protoplanet eccentricities than that of 
the disc. This will be a topic for future investigation. 

Although there is some uncertainty in the resonant 
torque calculations because of the need to sum contribu- 
tions of different sign, weakening of the tidal interaction is 
expected on general physical grounds on account of larger 
protoplanet disc relative velocities at higher eccentricity. 
This indication of migration reversal at the higher eccen- 
tricities was found to be supported by a local dynamical 
friction calculation applicable in that limit. In this case the 
interaction near apocentre tends to speed up the proto- 
planet while the interaction near pericentre tends to slow 
it down. These effects are of opposite sign but the longer 
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time spent near apoeentre results in a net outward migra- 
tion of the protoplanet for the surface density considered. 

Thus the existence of global non circular motions in 
discs with radial excursions comparable to or exceeding 
the 

semi-thickness may have important consequences for the 
migration of cores in the earth mass range. While pro- 
cesses of the type considered in this paper are unlikely 
to lead to the very high eccentricities observed for some 
giant planets, they may be important in controlling migra- 
tion during planet formation as well as producing modest 
eccentricities ~ 0.2. 

Acknowledgements. The author thanks the TAP for visitor sup- 
port and Caroline Terquem for valuable and stimulating discus- 
sions as well as a careful! reading of a preliminary draft of this 
paper. 
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The time averaged potential due to a perturbing 
inner planet 

The gravitational potential per unit mass due to a planet 
of mass nip located at rp = (r^, (pp) at r = (r, (f) is 



$ = - 



+r^ — ^rrp cos{9) 



(.64) 



with 6 = (fi — (fip. 

In order to incorporate the effects of protoplanets orbiting 
interior to the disc we adopt a Jacobi coordinate system. 
In this system the coordinates of the innermost proto- 
planet are referred to the central star. The coordinates of 
the remainder are referred to the centre of gravity of the 
central mass and all interior protoplanets. The disc is re- 
ferred to the centre of mass of central star and all inner 
protoplanets. 

This has the following consequences: 
For an object interior to the protoplanet with r < rp, 
the acceleration of the coordinate system due to the pro- 
toplanet must be allowed for. This gives rise, correct to 
first order in rUp, to the indirect potential ( eg. Brouwer 
& Clemence 1961) 



Grupr cos{9) 



(.65) 



to be added to the potential <&. 

For an object exterior to the protoplanet one must take ac- 
count of the fact that the coordinate system is now based 
on the centre of mass of the inner protoplanets and central 
star. Assuming initially that nip is the only such proto- 
planet, the central potential is modified to become 



$ = -. 



r -I- mpTp/M^l 



To first order in nip this gives 



GM* GnipTp cos 9 



(.66) 



(.67) 



In this case too the additional potential on the right hand 
side of the above can be incorporated into the perturbing 
planet potential and viewed as giving rise to an indirect 
potential. Thus the form of the perturbing potential 



$ = - 



Gnir, 



GmpTrp cos 9 



^r^+rl-2rrpCos{9) rnax{r^,rl,) ' 



(.68) 



APPENDIX 1 



incorporates both cases r < Vp amd Vp < r. 
In addition, although we included just one inner perturb- 
ing protoplanet, because we work to first order in their 
masses, the principle of linear superposition is valid such 
that the contributions of many such objects may be lin- 
early superposed. 

The perturbing potential due to a single protoplanet 
may be decomposed as a Fourier expansion in 6. 
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Thus 



m— 



COS m9. 



with 



1 



Grrip 7o 



27T 



<i>cos mOdO. 



(.69) 



(.70) 



From equation ( /75 ) the time averaged potential due to a 
protoplanet in circular orbit (^r(^p) = 0) is 



Ko{r,rp). 



(.76) 



When (rirp) ^ 0, equation ( |.75| ) gives the m — 1 com- 
ponent of the perturbing potential as the real part of 
$p(r) exp(i(^), with 



For the problem on hand, namely the study of global 

m = 1 modes of the disc planet system, we need only ^'„ir) = ;r— r 

consider m = and m = 1. Then from the above p 



Gnip^rirp) d{rlKi{r,rp)) 



drr. 



(.77) 



7r<i> 1 



Gnip 2 
with 



Ko{r,rp) 



Koir,rp) - Ki{r,rp)cos9, 



y^r^ + — 2rrp cos{6) 



de 



(.71) 



(.72) 



APPENDIX 2 



The equations of motion in two dimensional or- 
thogonal cordinates 

The basic equations of motion (|l]) written in vector form 
are 



and 



<9v ^ 1 

- + v.vv.--vn 



(.78) 



Kiir,rp) 




y^7'2 + r2 — 2rrp cos{9) 



cost 



■max{r^ , r-^) 



(.73) 



Time averaged potential for a protoplanet with 
small eccentricity 

We now suppose the protoplanet has a small eccentricity, 
Bp and write for its motion rp = ap{l + Cp cos tot), ipp = 
tot — 2ep sin ujt, where without loss of generality we assume 
the apsidal line to be along (p — where the displacement 
^r(^p) ~ Bp^p at ^ = 0. 

Here the protoplanet semi-major axis is Op and the orbital 
frequency is to. 

Expanding to first order in Cp, we find for the single pro- 
toplanet perturbing potential 



7r<I> 1 apCp dKo{r, ap) 
^ = --i^o(.,ap)-^ Q^cosut- 

dKi(r,ap) \ 
Ki[r,ap) + UpCp —cosbjt 1 cos((^— a;t-f-2ep sincji), 

CfCLp J 

Performing the time average then gives 

7r^o{r,ap) - ^^f- cosip. (.74) 



We wish to write these in component form in a gen- 
eral two dimensional orthogonal coordinate system (a. A) 
in the Cartesian {x, y) plane as introduced in section (||). 
The unit vectors in these orthogonal coordinate directions 
are ia = Va/|Va| and iA = VA/|VA| respectively. The 
orthogonal vertical coordinate z has an additional associ- 
ated orthogonal unit vector k = i^ x iA. However, there is 
no dependence on z and the velocity component in that 
direction can be taken to be zero. 
Thus we may write v = (uq, v\) or 



V = Vnln 



(.79) 



To deal with equation (.78), we first use the vector identity 



V- Vv- -V(|v| 



- a; X V, 



(.80) 



where a; = V x v is the vorticity (Arfken & Weber 2000). 

Because we only wish to consider the equations in a 
two dimensional limit , only the z component o f vo rticity, 
Lo^ is non zero so that we may write equation ([T^) as 



^ -f c^.k X V + -v(|vn = --VH - va>. 



(.81) 



We may now write equation (.81 ) in component form 
using the identities (see Arfken & Weber 2000) 



Grap 2n 



2apTT dap 



V = Va— +VA— 
oa OA 



For small Bp we may replace Op by rp and expressing the 
result in terms of the radial displacement ^ri^p), we find 



Grrir, 



Ko{r,rp) - i^^- cos (p. (.75) 



27r 



2r27r 



and 



(V X v) • k = |VA||Va| 



da 



dX 



(.82) 



(.83) 
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Doing this we obtain 



and 

^+c....+-|VA| 



1 ,<9n ,(9$ , 



After inserting (.83) into (.84) and (.85) it is a sim- 
ple matter to obtain equations (|34D and (|35|) as given in 
section (||). 



